#!/usr/bin/env python3
# coding: utf-8
"""
This module contains classes that represents a complete force fields and provides
methods in order to export these forcefields in various format as input for
molecular dynamics programs.
* ``ForceField`` is a general base class
* the ``ClassicalForceField`` class represents a classical force field defined from
a list of bonds, angles, dihedrals, impropers and non-bonded terms.
"""
__all__ = ['ForceField', 'ClassicalForceField']
[docs]class ForceField:
""" A general class to store a Forcefield """
def __init__(self, molecule):
"""
Args:
molecule (Molecule): The molecule associated to the force ForceField
"""
self.molecule = molecule
[docs]class ClassicalForceField(ForceField):
"""
A class which represents a classical force field such as Charm, AMBER, or
Gromos.
"""
def __init__(self, molecule, bonds=None, angles=None, dihedrals=None,
impropers=None, charges=None):
"""
All data needed to build the topology of the forcefield are mandatory.
Default values are empty lists.
Args:
molecule (Molecule): the Molecule
bonds (Bond): List of bonds
angles (Angle): List of angles
dihedrals (Dihedrals): List of dihedral angles
impropers (Improper): List of improper torsions
charges (list): atomic charges
"""
super().__init__(molecule=molecule)
self.bonds = bonds
self.angles = angles
self.dihedrals = dihedrals
self.impropers = impropers
self.charges = charges
# TODO: on pourrait imaginer une méhtode qui retourne l'énergie ???
[docs] def as_dict(self, params="all"):
"""
Return a serializable dict to export whole or part of the force field
to a json file. Look at monty.json and MSNoable class for the future.
The approach is the one followed by pymatgen.
Args:
params (list or string): List of forcefield component to be exported
among 'bonds', 'angles', 'dihedrals', 'impropers', 'charges',
'vdw'. If 'all', all the forcefield is return.
"""
if params == "all":
export = ["bonds", "angles", "dihedrals", "impropers", "charges", "vdw"]
else:
export = params
forcefield = {"molecule": self.molecule.as_dict()}
if "bonds" in export:
forcefield["bonds"] = [bond.as_dict() for bond in self.bonds]
if "angles" in export:
forcefield["angles"] = [ang.as_dict() for ang in self.angles]
if "dihedrals" in export:
forcefield["dihedrals"] = [dihe.as_dict() for dihe in self.dihedrals]
if "impropers" in export:
forcefield["impropers"] = [imp.as_dict() for imp in self.impropers]
if "charges" in export:
forcefield["charges"] = [q for q in self.charges] if self.charges else None
if "vdw" in export:
# Not supported yet
pass
return forcefield
[docs] def export_lammps_data(self, filename=None):
""" Export the force field for LAMMPS """
def write_parabolic_coeffs(coords, title):
""" write coefficients in case of parabolic potential """
lines = ""
if len(coords) > 0:
lines += "\n %s\n\n" % title
for ic, coord in enumerate(coords, 1):
lines += " %2d %12.1f %6.2f\n" % (ic, coord.frc_cste, coord.value)
return lines
def write_internal_coordinates(coords, title):
""" Write the list of internal coordinates coords """
lines = "\n %s\n\n" % title
for ic, coord in enumerate(coords, 1):
lines += " %2d %2d " % (ic, ic)
lines += " ".join("%3d" % (iat + 1) for iat in coord.atoms)
lines += " # "
lines += " ".join("%2s" % el.specie.symbol for el in coord.elements)
lines += "\n"
return lines
lines = 'LAMMPS data file generated by MAMMOTH version 0.1\n'
lines += '%d atoms\n' % len(self.molecule)
lines += '%d bonds\n' % len(self.bonds)
lines += '%d angles\n' % len(self.angles)
lines += '%d dihedrals\n' % len(self.dihedrals)
lines += '%d impropers\n' % len(self.impropers)
lines += '%d atom types\n' % len(self.molecule)
lines += '%d bond types\n' % len(self.bonds)
lines += '%d angle types\n' % len(self.angles)
lines += '%d dihedral types\n' % len(self.dihedrals)
lines += '%d improper types\n' % len(self.impropers)
# TODO: il faut vraiment un tab ?
cart_coords = self.molecule.cart_coords
lines += '%6.2f %6.2f \t xlo xhi\n' % (3 * cart_coords[:, 0].min(),
3 * cart_coords[:, 0].max())
lines += '%6.2f %6.2f \t ylo yhi\n' % (3 * cart_coords[:, 1].min(),
3 * cart_coords[:, 1].max())
lines += '%6.2f %6.2f \t zlo zhi\n' % (3 * cart_coords[:, 2].min(),
3 * cart_coords[:, 2].max())
lines += write_parabolic_coeffs(self.bonds, "Bond coeffs")
lines += write_parabolic_coeffs(self.angles, "Angle coeffs")
if len(self.dihedrals) > 0:
lines += '\n Dihedral Coeffs\n\n'
for idih, dihedral in enumerate(self.dihedrals, 1):
if 0 < abs(dihedral.value) <= 90:
phi0 = 0
elif 90 < abs(dihedral.value) <= 270:
phi0 = 180.
elif 270 < abs(dihedral.value) <= 360:
phi0 = 0
# TODO vérifier si cette convention de phase est bonne
# sachant que 0 est cis et 180 est trans
lines += " %2d %8.1f %2d %6.1f %6.1f\n" % (idih, dihedral.frc_cste,
dihedral.periodicity,
phi0, 1.0)
lines += write_parabolic_coeffs(self.impropers, "Improper coeffs")
lines += '\n Masses\n\n'
for a, sp in enumerate(self.molecule.species, 1):
lines += '%3d\t%7.3f # %s\n' % (a, sp.atomic_mass, sp.symbol)
lines += '\n Atoms \n\n'
for a, site in enumerate(self.molecule):
data = (a + 1, a + 1, self.molecule.at_charges[a],
site.coords[0], site.coords[1], site.coords[2], site.specie)
lines += "%3d 1 %3s %9.3f %10.4f %10.4f %10.4f # %2s\n" % data
lines += write_internal_coordinates(self.bonds, "Bonds")
lines += write_internal_coordinates(self.angles, "Angles")
lines += write_internal_coordinates(self.dihedrals, "Dihedrals")
lines += write_internal_coordinates(self.impropers, "Impropers")
lines += "\n"
if filename:
with open(filename, "w", encoding="utf8") as f:
f.write(lines)
else:
return lines
def __str__(self):
""" Return a str view of all the data """
lines = "FORCE FIELD DATA\n"
lines += "================\n\n"
lines += "Bonds\n"
lines += "-----\n"
for bond in self.bonds:
data = (bond.atoms + tuple(el.specie for el in bond.elements) +
(bond.value, bond.frc_cste))
lines += "%4d %4d (%2s, %2s) %18.3f %8.1f\n" % data
if self.angles:
lines += "\n"
lines += "Angles\n"
lines += "------\n"
for angle in self.angles:
data = (angle.atoms + tuple(el.specie for el in angle.elements) +
(angle.value, angle.frc_cste))
lines += "%4d %4d %4d (%2s, %2s, %2s) %14.3f %8.1f\n" % data
if self.dihedrals:
lines += "\n"
lines += "Dihedrals\n"
lines += "---------\n"
for dihedral in self.dihedrals:
data = (dihedral.atoms + tuple(el.specie for el in dihedral.elements) +
(dihedral.value, dihedral.frc_cste))
lines += "%4d %4d %4d %4d (%2s, %2s, %2s, %2s) %10.3f %8.1f\n" % data
if self.impropers:
lines += "\n"
lines += "Impropers\n"
lines += "---------\n"
for improper in self.impropers:
data = (improper.atoms + tuple(el.specie for el in improper.elements) +
(improper.value, improper.frc_cste))
lines += "%4d %4d %4d %4d (%2s, %2s, %2s, %2s) %8.3f %8.1f\n" % data
return lines